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data is analyzed. The method developed in this Paper combines the power 
of the Proper Orthogonal Decomposition (POD) technique to identify these 
coherent structures in experimental data sets, and its optimality in providing 
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the procedure to analyze coherent structures in an oceanic setting, the ones 



POD basis used is the one obtained from the experimental data. We apply 



arising from instabilities of the Algerian current, in the western Mediterranean 
Sea. Data are from satellite altimetry providing Sea Surface Height, and the 



^ ■ model is a two-layer quasigeostrophic system. A four-dimensional dynamical 

system is obtained that correctly describe the observed coherent structures 
(moving eddies). Finally, a bifurcation analysis is performed on the reduced 
model. 
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I. INTRODUCTION 

In the last decades the study of turbulent or extended chaotic systems has enjoyed 
important advances. Two of them are, first, the recognition of the existence and high 
relevance of coherent structures (defined as strongly persistent spatiotemporal structures 
in the dynamical evolution of the system) in weakly and even strongly chaotic systems, 
and, second, the borrowing of mathematical methods coming from the studies of nonlinear 
dynamical systems (see PJ and references therein). 

In both subjects, the introduction of the statistical technique known as the Proper 
Orthogonal Decomposition (POD, also known under a variety of other names, such as 
Karhunen-Loeve decomposition, method of Empirical Orthogonal Eigenf unctions, etc.) has 
played an important role. It was introduced in the context of turbulence by Lumley || and 
has revealed itself as an efficient technique for finding, describing and modeling coherent 
structures in turbulent fluids or extended chaotic systems. The purpose of POD is to sep- 
arate a given data set into orthogonal spatial and temporal modes which most efficiently 
absorb the variability of the data set. 

The power of POD has been implemented following two different paths QXJ] : On the 
one hand the POD is used as a standard technique to extract coherent structures from 
empirical data sets @,[§]. Contrasting to Fourier Decomposition, the eigenfunctions obtained 
from the POD may display spatial localization, and thus provide a more efficient way to 
represent coherent structures. The use of empirical information can be pushed further and 
methodologies from dynamical systems theory and other fields have been used in the POD 
framework, to provide useful algorithms for control || and prediction |||7|. On the other 
hand, the POD eigenfunctions provide a set of basis functions which is optimum (at least 
in a well defined linear sense) for obtaining low-dimensional ordinary differential equation 
(ODE) approximations starting from models based on partial differential equations (PDEs). 
The approximation is performed by obtaining long runs of the PDEs, performing the POD 
onto this synthetic data set, and using the Galerkin method to project the PDE model into 



the so obtained POD eigenfunctions fq-fl2|l- 

These two potentialities, i.e. the ability to extract empirical information from experi- 
mental data, and the efficiency in building low-dimensional projections from models, are not 
frequently used together in the literature. A remarkable exception is the use of empirical 
eigenfunctions obtained from the POD of experimental data from a turbulent boundary 
layer to build a low- dimensional approximation to the Navier-Stokes equations |13| , [T^| . We 
believe that, in some circumstances, the projection of theoretical models into experimentally 
obtained empirical functions could improve both the model and the data. This will occur 
in situations such as in the modeling of natural phenomena (ocean or atmospheric dynam- 
ics, for example) where even very complex models may be not accurate enough, and data 
are unvoidably noisy and difficult to calibrate. Projecting the model onto the experimental 
eigenfunctions will force it to stay into the 'right' subspace, providing a kind of data assimi- 



lation [15| that may compensate the loss of details inherent to low-dimensional projections. 
On the other hand, the truncation involved in the POD method implies a kind of filtering 
providing noise reduction to the data set. 

Our aim in this Paper is to explore the synergy between experimental observation and 
low-dimensional reduction via POD, in the complex setting of environmental fluid dynamics. 
In particular, a model for coherent structures arising from instabilities of the Algerian current 
in the Mediterranean Sea will be set up and analyzed. 

In the ocean dynamics context, the kind of spatiotemporal data sets we need for our pur- 
poses can only be obtained from satellite observations. The recent availability of satellite 
data of the sea surface is allowing a deeper understanding of the Ocean. Satellites contin- 
uously measure sea temperature, sea level, chlorophyll concentration, etc., which increase 
our knowledge of ocean currents, mean sea level changes, tides, or plankton dynamics, to 
name a few. In particular, in the last decade the ERS and the TOPEX/POSEIDON (T/P) 
satellite missions have provided the scientific community with high-accuracy altimetry data, 
which determines the sea surface level, this is, the height of the sea surface over a refer- 
ence level on Earth. Among other relevant scientific applications, these type of data are 



specially useful for a better understanding of the dynamics of mesoscale phenomena in the 
Ocean. Mesoscale refers to typical spatial scales of 30 to 300 kilometers and time scales of 
less than one year, and it is associated with movements of oceanic currents and short-time 
flow variations, and also with the formation and propagation of ocean eddies. These eddies 
are generated by interactions with the oceanic topography and/or mean flow instabilities 
and their importance is enormous: for example they play a fundamental role in the heat 
transport from low to high latitudes. 

Some details of the results we present here are determined by the peculiarities of the 
data set we are going to use. In particular we mainly focus in the motion of a vortex which 
is present in the data, but is only described by subdominant eigenf unctions in the POD 
results. This lead us to a model for this coherent structure, but we do not try to model the 
full dynamics of whole data set. We expect however that our general methodology will be 
useful in other problems in which real-world noisy observations and complex but imperfect 
PDE models are available. Generalizations of the POD method which take into account in 
a more consistent way dynamic constraints have been developed flT5|-|r8[ and even applied to 
geophysical contexts |TB],[TP],|2D| . We will use however, and just for simplicity, the standard 
formulation of the POD technique. 

The Paper is organized as follows: in the next section, we describe the data set. In section 
III, these data are analyzed with the Proper Orthogonal Decomposition (POD). The study 
of the most relevant of the eigenfunctions allows the identification of a coherent structure, 
that is, a moving vortex or eddy. Then, in the following section, the associated temporal 
modes of the POD eigenfunctions defining the eddy are projected over a hydrodynamic model 
which, finally, provides a deterministic dynamical system depending on the parameters of 
the model. In section V, the reconstruction of the moving vortex from the dynamical system 
model is performed. Next, in section VI the bifurcation analysis of the dynamical system is 
shown. Section VII concludes this work. 



II. SATELLITE ALTIMETRY DATA 



We analyze altimetry data from the T/P and ERS-1 satellite missions [21]. Altimetry 
data of the Ocean provide the Sea Surface Height (SSH) over a reference substrate. The data 
obtained in both missions have been merged (to obtain a better spatiotemporal resolution) 
on a common time period, from October 1992 to December 1993 and, finally, 44 maps taken 
every 10 days on a 0.2° regular grid are obtained for the Western Mediterranean Sea [^5|. We 
restrict our analysis to the area known as the Algerian Current, localized between — 15° E 
and 35 — 40°iV, where a strong mesoscale activity is observed [24|. A mean flow moving 



eastwards and parallel to the coast of Algeria is the main feature in this area. It undergoes 
instabilities that shed vortices into the western Mediterranean basin, greatly influencing the 



physical and biological processes in this area of the Sea [23| . In Figure 1 a) we show in a 
small box of an image of the Mediterranean Sea, the area under study. Figure 1 b) shows 
one of the 44 maps that we are going to analyse. 

The SSH fields allow the identification of coherent structures (structures approximately 
maintained in the flow over long time periods) in geophysical flows. In particular, areas 
of higher altimetric values may correspond to anticyclonic (clock-wise) vortices and lower 
ones may indicate the existence of cyclonic (anticlock-wise) vortices. Actually, the data we 
have used are referred to a mean level, i.e, we analyze Sea Level Anomalies (SLA) where 
the reference height is the temporal mean of the data. This may give rise to some minor 
problems because to obtaining the SSH data (the one we are going to need in our modeling 
approach) is not as simple as adding the mean sea level. This is because of the different 
resolution in the data and will be explained in detail in section IV. 

III. POD ANALYSIS OF THE SATELLITE DATA 

As it has already been mentioned, the POD technique is generally used to analyze ex- 
perimental or numerical data with the view in extracting their dominant features, which 



will typically be patterns in space and time. On output, it provides a set of orthogonal 
functions which are the eigenfunctions of the covariance matrix of the data. Generally, this 
set is ordered in decreasing size of the corresponding eigenvalue, the larger the eigenvalue 
meaning the larger percentage of the data variance is contained in the dynamics of corre- 
sponding eigenfunction. Thus, if u(x,y,t) is our data field ((x,y) G A C R 2 is the spatial 
point and t is time) to which the temporal average has been subtracted, the POD basis 
{4>i(x,y),i = 1, ...,00} is obtained after solving 



< u(x, y, t)u(x', y', t) > 4>i(x', y')dx'dy' = Xi4>i(x, y), (1) 

I A 

being < • >= i / • dt, i.e., the time average, and A, the corresponding eigenvalues, which 

are ordered in decreasing size Ai > A2 > ... > 0. Therefore, we have the modal decomposition 

00 
u(x, y,t) = J2 ai{t)4>i{x, y), (2) 

i=i 

where the a«(t) are the so-called temporal modes. In addition, the following orthogonality 

conditions are fulfilled 

<j>k(x,y)<j>i(x,y)dxdy = Ski, (3) 

J A 

1 f T 

(a k (t)ai(t)) = — / a k {t)ai(t)dt = \ k 5 kh (4) 

1 Jo 

where 8 k i is the Kronecker delta. 

The optimality of the POD basis functions means that Jl[, among all linear decom- 
positions with respect to an arbitrary basis {3>i(x, y)}, for a truncation or order N, i.e., 
u N (x,y,t) = ^2iLiai(t)^i(x,y), with a,i(t) = f A u(x,y,t)&i(x,y), the minimum error, defin- 
ing the error as e = (j A (u — u N ) 2 dxdy\ is obtained when {&i(x, y)}is the POD basis 
{(j)i(x,y)}. 

We now apply the POD analysis to the altimetry satellite data. The results of this 
are outlined in the following. Fig. |^ shows (in linear-log scale) the fraction of variance 
Aj/(Z)^=i A„), given by each eigenvalue. It is clearly seen that most of the variance is 
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captured by the first and second eigenvalues. In Fig. [|we show the temporal mode associated 
to the first two eigenvalues, and the power spectrum of the four dominant ones is plotted 
in Fig. |j. The annual seasonal periodicity is clearly observed in the two dominant modes. 
This and more detailed observations in terms of a complex version of the POD on the same 
data set p3] allow us to interpret the dynamics given by these dominant eigenfuntions as 



the seasonal response of the Ocean (heating in summer and cooling in winter) along the 
annual cycle. 

Therefore, to get a deeper insight into the data (and in particular in the mesoscale 
phenomena), we have to study other eigenfunctions than the first two. In Fig. |5] we show 
again the fraction of variance of the eigenvalues (linear-log plot) but in this case we have 
removed the first two eigenvalues. We observe that eigenvalues 3 rd and 4 th are equivalent 
under the error bar (for calculating error bars in POD eigenvalues see ||25|| ). In this case, 
a linear superposition of them can represent a moving coherent structure ||. In Fig. |6| 
we show the temporal modes associated with eigenvalues 3 rd and 4 th . This Figure, and the 
corresponding power spectra in Fig. ^ suggest a weak semiannual periodicity for both modes. 

Visualization of the time evolution of the data filtered to keep just the eigenfuntions 3 rd 
and 4 th , i.e., uu(x,y,t) = a 3 (t)0 3 (x, y) + a^(t)(p4(x,y), suggest a vortex moving northward 



and eastwards from the Algerian coast, in agreement with the results of ||24| . The approxi- 
mate periodicity of the associated temporal modes indicates that new vortices are shed from 
the coast roughly every six months. Describing the dynamics of this coherent structure 
found in the flow will be our goal in the remaining of the Paper. 

IV. MODEL AND LOW-DIMENSIONAL DYNAMICAL SYSTEM 

The classical way of using POD to obtain a low-dimensional dynamical system approxi- 
mation, which takes into account the most relevant features of the physical system, consists 
in truncating the expansion (|2[) to a particular order |T|,|5|-|H|] . This order is generally chosen 



to contain most of the percentage of the variance of the data. Then, the equations gov- 



erning the dynamics of the system (PDEs from which the data may have been generated 
numerically) are projected over this particular Galerkin basis and a system of ODEs for the 
temporal modes di(t) can be obtained. Our approach is somewhat different, first of all, our 
data are from satellite observations and we need a specific mathematical model to describe 
approximately our data, and second, our interest focuses in the dynamics associated with 
the eigenfuntions 3 rd and 4 th which seem to contain the evolution of the moving mesoscale 
vortex, and not the more dominant eigenfunctions 1 st and 2 nd . 

Proceeding with the modeling step, and supported by marine experimental campaigns 



24] , |26|1 , we assume that the strong mesoscale activity in the Algerian current is mainly due 
to baroclinic instability phenomena. This name refers to the instabilities grown from the 
available potential energy associated with horizontal gradients of density |[28| . Therefore, 



we choose a two-layer quasigeostrophic model as our basic flow description since this is the 
minimal model accounting for these type of instabilities. Layered quasigeostrophic models 
are widely used in oceanographical modeling and their main assumption is that the Ocean 
behaves as having different layers where density is constant and, in all the different layers, 
geostrophic balance is maintained (i.e. Coriolis and pressure forces nearly equilibrate via a 
quasibidimensional flow). Actually, the spatial scales of the chosen region, and its strong 
topographic features, lead however to important deviations from quasigeostrophy. Thus, 
the postulated model should be considered at most as a crude approximation to the real 
dynamics. It is one of the objectives of this Paper to show that the the empirical information 
contained in the satellite data is incorporated into the model during the projection procedure, 
so that the final low- dimensional model gives a reasonable description of the dynamics. 

In the framework of multilayer quasigeostrophic models, every fluid layer i of density pi 
and thickness hi is described by a stream function tpi, which is proportional to the pressure 
field within the layer, and such that the horizontal velocities Vi = (ui, Vi) within the layer 
verify Ui = — q and Vi = -4^. In our equations, the coordinate directions x and y will be 
oriented along the northward and the eastward directions, respectively. 

More specifically we use a two-layer quasigeostrophic model on a beta plane and over 
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topography. An eddy-viscosity and a bottom friction terms are also included. The equations 
defining the dynamics of the stream function of both layers are 

V> 2 - Vl 



Ei 

Dt 



V^i + / + 



m 



Ei 

Dt 



V ^2 + / + — ^ — + /ott 

-"2 



i? 2 . 



z/v> 2 - c 6 v 2 v> 2 , 



(5) 
(6) 



where the subscript i = 1 (2) refers to the upper (bottom) layer, -^ = -^ + J(if) i: •), ijji(x,y,t) 
is the layer stream function, r(x, y) is the bottom topography, Ri 



M^,R, 



N 2 H 2 
/o ; 



Nf = -§-^-, N 2 = ^- e -, where Sp — p 2 — p\. g the gravitational acceleration, Hi is the 
mean thickness of the layer i, and / is the Coriolis parameter, which in the beta-plane 
approximation depends on the latitude as / = f + (3y. v is the eddy-viscosity and Cb is the 
coefficient describing friction with the bottom of the sea. The Jacobian operator J(A, B) is 
defined as: 



_dAdB__dAdB_ 
dx dy dy dx 



(7) 



More details about quasigeostrophic dynamics can found for example in Refs. [^] and 
. A more explicit way to write our equations (0) and 



is 



dV 2 ^ , 1 d(i/j 2 
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+ P 



dj> 2 
dx 



(9) 



In this model the stream function of the upper layer is tpi(x,y,t) = j-h(x,y,t) where 
h(x, y, t) is the height of the sea surface over the point (x, y) at time t. This last quantity 
is the one linked to the satellite observations on which we have performed the POD. It is 
important to note that in equations (§) and (®) we have not considered an annual forcing 
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term which accounts for the the seasonal variability of the Algerian Current. This is a very 
important fact for the next step in our approach, the projection onto a particular Galerkin 
basis determined from the observations. We assume the following ansatz: 

if}i(x,y,t) = < ip(x,y,t) > +$'(x,y,t) 

= — (< h(x, y, t) > +a 3 (t)4) 3 (x, y) + o 4 (f)<M^ y)), (10) 

Jo 

where the temporal coefficients 03(f) and 04(f) will be calculated in the following sections and, 
finally, compared with the POD temporal coefficients 03(f) and 04(f). With this we assume 
that the stream function %p\ (which is proportional to the height) of the upper layer is de- 
composed in its temporal mean < ip(x,y,t) >= 4- < h(x,y,t) >, accounting for the annual 
mean flow, and a perturbation ip'(x,y,t) = j- (03(f) 03 (a;, y) + 04(^4(0;, y)), which models 
the mesoscale processes [^7fl. The perturbation basis for the height {03(x, y), 4>4(x,y)} is 
what we have obtained from the POD analysis of the data, and the way to calculate the 
annual mean will be detailed at the end of this Section. 

As we have no real measure for the bottom layer (the satellite sensors get data just from 
the sea surface), we need to make some additional hypothesis in our model. We propose the 
following ansatz for the expansion of the bottom layer's stream function 

i> 2 {x, y, f) = -Uy + Vx + |-(6 3 (f)0 3 (x, y) + b 4 (t)fa(x, y)), (11) 

Jo 

where U and V are parameters of our model and simulate the eastward and northward 
velocity, respectively, of the bottom layer flow. The physical meaning of ( |TTD is that the 
perturbations from the mean flow for the bottom layer are generated by the same basis 
functions as the upper one though with, obviously, different temporal coefficients. The most 
important feature in the former ansatz is the mean flow, parameterized with U and V. The 
values of these parameters not only determine the intensity of the mean flow but also, and 
most importantly, its direction and sense. Discussions about the role of the different values 
of U and V will be given in the next section. 

Projecting expansions (|T0D and (|TT|) over equations (g) and ([]) and using the orthogo- 
nality relations of the POD basis 4>%{x, y), we obtain the evolution equations for the coherent 
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structure's temporal amplitudes a, and b% (i = 3,4), 

+ P6,i b 3 + P7,i b l + PS,ih + P9,ih + Pio,ihh + Pll,i, 

-77 = qi,ia>3 + Q2,ial + g 3 ,j03 + Q4,id 2 + q5,ia 3 a 4 

+ qe,i b l + <l7,i b l + qs,ih + q9,ih + qw,ihh + qn,i- (12) 

The coefficients pu,i and qk,i (k = 1, .., 11 and i = 3,4) are real numbers that depend on 
the parameters of the model and on integrals containing < ip >, 3 (x, y), 4>^(x,y), the and 
their derivatives. Their explicit expressions are quite complex and have been obtained by 
computer algebraic manipulation. We do not write down here all these involved expressions. 
Just to give an example of them, in Appendix A we display the mathematical expression for 

9l0,2- 

We finally proceed to explain how to obtain the mean flow < ip(x, y,t) >, without which 
the coefficients in ( |T2] ) remain undetermined. Some manipulations of the data are needed 
because of the bad spatial resolution of the available mean field. We recall that we are dealing 
with Sea Level Anomaly (SLA) data obtained from the two altimetric missions ERS-1 and 
T/P. These are conveniently treated to obtain regular maps in space and time every 10 
days and on a 0.2° regular grid. These SLA are relative to the annual mean sea level and, 
therefore, we need this annual mean to obtain the total height of the sea surface and thus be 
able to make relation between the empirical eigenf unctions 03, 04 and the dynamic variable 
•01 of the quasigeostrophic model (see Eq. (T0|) . Unfortunately, the annual mean sea level 
data are not manipulated to improve their resolution, and we have just the T/P data, of a 
very coarse resolution (around 2.8°) to calculate this annual mean. Therefore, we need to 
interpolate these to a 0.2° regular grid and then to add the resulting annual mean to the 
SLA data, in order to obtain a consistent SSH field. But, all these manipulations are, at the 
end, manifesting when we solve eq.(|T2|) in such a way that a 3 (t) and 04 (t) have a nonzero 
average, in contrast with the POD temporal modes, a 3 (t) and a^t), calculated from data, 
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- / Oi{t)dt^Q } 1 = 3,4. (13) 

1 Jo 

In order to heal this, we proceed with an assimilation-like approach: thus we modify the 

annual mean sea level, which has been obtained interpolating the T/P data, by adding this 

nonzero average, i.e., 

f 1 r T 
h m {x,y) = h T [ p (x,y) + — — / a 3 (t)((j) 3 (x,y,t) + a 4 (t)(fr 4 (x,y,t))dt, (14) 

g 1 Jo 

being h m (x,y) the new annual mean height and h^/ p (x,y) is the interpolated T/P annual 

mean height. Finally, and with high numerical accuracy, the new temporal modes obtained 

from eq. fll2"|) have now temporal zero average. 

Summing up, the data we are using along this paper are obtained by adding to the SLA 
data the above calculated h m (x,y) field. It is important to note that because in the POD 
analysis we substract the mean field of the data, the qualitative features of the analysis 
are similar if we analyse the SLA data or the SLA plus the h m field. Nonetheless, the 
quantitative differences are not negligible at all at the level of the dynamical system (|i~2"l). 

The four-dimensional dynamical system (|i~2"D, now fully defined, is the desired low- 
dimensional approximation aimed to describe the coherent structures in our data set. In the 
next Section we show that the dynamics of the observed coherent eddy is recovered from 
(0). 

V. NUMERICALLY GENERATED COHERENT STRUCTURE DYNAMICS 

We now proceed to integrate the equations (0). First, typical values for the parameters 
of the quasigeostrophic model (|]) and (Jffy are needed. At mid-latitudes, adequate values for 
the parameters giving the Coriolis force are (3 = 1.0 10~ n m/s and fo = 10~ 4 s _1 f28| , |3"0"|j . In 
addition, for the Algerian Current area the values p\ = 1025 kg/m 3 and pi = 1029 kg/m? 
are adequate for the densities of the upper and bottom layer and H\ = 300 m, H2 = 3500 m 
for their mean heights. The bottom topography t(x, y) are real data obtained from the data 
basis at the URL ' |http://modb. oce.ulg.ac.be/Bathymetryhtmlf . For the friction parameter 
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with the bottom topography we take a value of C& = 1.5 10~ 6 . A typical value for the eddy 
viscosity at the scales we are working is v — 200 m 2 /s. The remaining parameters are the 
geostrophic velocities of the bottom layer. This is a difficult task since nobody really knows 
what happens in the deep waters of the Algerian Current. Anyway, recent results obtained 
by the PRIMO-1 experiment in the channel of Sardinia |23| seem to indicate that typical 
velocities for the deep layer are rather weak (of the order of 1 cm/ s), though it is not clear 
if the deep layer flow (dlf ) is westward or eastward. Therefore we assume dlf following the 
upper layer, this is northeastward, with typical values of: U = 5 cm/ s and V = 5 cm/ s. In 
the next Section, other values of U and V will be discussed. 

Fig. [7] shows a 3 (t) and a 4 (t) obtained by integrating ( |l"2"D with a fourth order Runge- 
Kutta method. The periodicity of both is evident, being the period of around 5 moths, 
in good agreement with the observed one. If some annual forcing would be added to our 
model, this period would probably lock to the semiannual harmonic, still improving the 
agreement. The shape of the oscillation in the calculated evolution is much more regular 
than the experimental one, as expected from a clean simulation versus noisy data. In Fig. |8| 
we show a temporal sequence of the coherent structure dynamics given by a^{t)(j)^{x,y) + 
a4(t)0 4 (x, y) from the calculated temporal modes. An eddy appears next to the coast and 
moves northeastward, the process being repeated five months later. This is fully consistent 
with the observed data and confirms the success in our task of obtaining a low-dimensional 
reduced model. Sensibility of the results to variations in the values chosen for the parameters 
is discussed in terms of a bifurcation analysis in the next Section. 

VI. BIFURCATION ANALYSIS 

The value of the eddy viscosity v is somehow arbitrary since it would depend on the scale 
of observation. Thus, a discussion of the variations in model behavior as v varies is in order. 
We analyze numerically our system of four ODE's (|1^) with the help of the software package 
Dstool [[TT| which integrates the system with a fourth order Runge-Kutta. We change the 
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eddy-viscosity parameter v, which is proportional to the inverse of the Reynolds number, 
and observe the bifurcation behavior. The rest of the parameters take the same values as in 
the former section. It should be noted that in principle the empirical eigenfunctions would 
vary in a system with varying u, but since we just have experimental data for the actual 
value of the eddy viscosity in the Ocean at the observed scales, we keep the parameters in 
(|l2l) as determined from the POD of the observed data. The bifurcation diagram is outlined 
as follows: 

• For v > 212 m 2 /s there are six fixed points. Two of them are stable and the rest are 
unstable. 

• When v = 212 m 2 /s a Hopf bifurcation occurs. One of the stable fixed points (the 
one localized at the origin) gets unstable by decreasing v and a limit cycle appears 
surrounding it. The limit cycle persists for all the values of the viscosity smaller than 
212 m 2 /s. The system undergoes no new bifurcations by decreasing the viscosity 
parameter. 

The main dynamical feature in our model is thus the existence of a Hopf bifurcation 
which gives birth to a limit cycle for a long range of v values, including the physical ones at 
the scales we are working. All these limit cycle solutions give rise, after reconstruction of the 
coherent structure with the help of the empirical eigenfunctions to traveling wave patters 
with a period of around six months. In particular, the moving eddy identified in section [V| 
is just one of these solutions. Moreover, the rest of the fixed points in the second regime, 
i.e. when v < 212 m 2 /s, seem to have no physical significance as their basins of attraction 
correspond to very high values of the initial condition for a.j(t) (i = 3, 4), i.e., high values of 
the sea surface height. When the eddy-viscosity is too large, the system evolves towards a 
stable fixed point, with no moving coherent structures, as expected on physical grounds. 

To give a stronger support to the former analysis, we have tested our ODE system with 
other values of the dlf velocity. We have observed numerically the following behaviour of 
the system: 
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• High values of U or V, ~ 10 cm/s, produce diverging solutions, with no fixed points 
for any value of v. Very low values (~ 10~ 2 cm/s) give rise also to unbounded solutions 
for any initial condition in some range of high viscosity values (low Reynolds number), 
which is not reasonable on physical grounds. 

• Changing sign in U, that is, assuming a westward flow in the bottom layer, does not 
change considerably the bifurcation diagram, but the amplitude of the limit cycle is 
too big when reasonable values of the viscosity (around 200 m 2 /s) are used. On the 
contrary, changing sign in V or in U and V simultaneously gives rise to a ODE system 
where all the solutions are unbounded. 

We think these are enough reasons supporting the chosen direction and magnitude of the 
velocity of the dlf, that is, a northeastward direction and typical values for the horizontal 
velocites around 5 cm/s. 

VII. SUMMARY 

In this article we have used the Proper Orthogonal Decomposition to obtain a low- 
dimensional dynamical description of coherent structures observed in satellite data of a 
region of the Mediterranean Sea. First, analysis of the altimetric satellite data via the POD 
allows the identification of a moving vortex in the ocean surface. Second, projection of 
a two-layer quasigeostrophic model onto the empirical basis, together with some physical 
assumptions on the unobserved part of the sea, allow the construction of a fourth-order 
dynamical system that gives a reasonable description of the dynamics of the coherent struc- 
ture, in particular its period and amplitude. It is remarkable that a crude PDE model, and 
noisy data, can be merged to obtain an efficient reduced model. We expect our general 
methodology would be of use in other complex environmental fluid dynamics problems. 
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VIII. APPENDIX 



If we denote the scalar product by {03, 04} = J A 03 (x, y)04(x, y)dxdy, and we define also: 
c k ,i = -{A(fi k ,(pi}, (15) 



, a d(f) k 4 g d <ip> dA<p k dA<ip>d<i) k 

> l = i~^ir~ + ^v fc - - — — + — — 

ox jo ox dy dy ox 

_ d < ip > dA<f) k _ dA<ip> d<p k 9_ (lJ d<h _ y ^\ ± 1 / lg \ 

dy dx dx dy foRf dx dy 

9 [ d<i])>d(j> k d<^>d<j) k 
f$R 2 dy dx dx dy ' 

fo f drd(l> k drd(f) k d$ k 4 

ii 2 c/y ox ax dy dx 

g , dA(f) k _ dA(f) k _ g d < ijj > d<p k _ d < ijj > d<p k 

/o dy dx / i?2 dy dx dx dy 

a . d&iu dchu 



r d(f) k d(f) k 
dv dx 



{vT~ u T,m, (is) 



foRi dy 

__ , g d <ip >dA <ip > g d < i/j > dA < ip > g <9 < -0 > 
/o dx dy f dy dx f R 2 dx 



(19) 



(20) 



T ,<9 < ^ > „d <ib > „ A 

+ Ai?' '-fi— l> -£- + '" 4 <*>•«• 

r g y g < ^ > ff 9 < -0 > J±tt<^ 

q2 ' 1 l f Rf dy + f Rf dx + H 2 U dx 

d g_| g03_gA03 90 3 9A03 

/o <9x <9y dy dx ' 

d g_. g04.gA04 d0 4 (9A0 4 

/o <9x cfy dy dx ' 

g a0 3 aA0 4 <90 3 <9A0 4 <90 4 <9A0 3 <90 4 <9A0 3 

^ = { "^ ( 9x-^T " "^^T + 5x"^T " ^^T' * ,} ' (23) 

with j, / = 3, 4, we finally obtain, 

910,2 = l - (-Rl -R2 ^3,3C4,4C3,4 — C3 3 G?3 4-R2 -Rl + -R2 <^3,3C3,4 + Ci o-Rl ^3,4-^2 C4,3 
fel V 

+ i?i i?2 £3,403,3 — C3 ; 3i?i (^3,4-^2 — Q,4-Rl C 3 , 3 d 3)4 R2 ~ C 4,4-Rl C 3 , 3 R2 ^3,4 

_ C 3,3-R2 <^3,4 + C34-R1 ^3,4-^2 C4,3C3,3 ~ £3,4 -Rl C4,3-R2 ^3,3 + C^Ri C 3y3 R2 C 3A d 3j3 

- c 4A Ri 4 c 3 /d 3A R2 4 + c 3 , 3 R 2 4 c 3A Ri 2 d 3i3 ) , (24) 



where 
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k\ = — 2c 4j3 -R 2 C 3j4 i?i + C 4>4 _R 2 R\ C 3j3 + C 4j4 _R 2 -Rl C 3> 3 + C 4)4 -R 2 -Rl C 3)3 

+ C 4j4 -Ri C 3j3 — C 4;4 i?i C 4[3 i?2 C 3)4 + 2c 4i4 /?i C 3i3 i?2 + C 4 ,4-Rl C 3j3 i?2 _ C 4)3 i?2 C 3i4 i?i C 3j3 
+ C 4 ,3 -R2 C 3;4 i?i — C 4j3 i?2 C 3i4 i?i C 3)3 — C 4i4 i?2 C 4 ,3 c 3,4-Rl + C 4 ,4-R2 C 3;3 + C 4i4 i?2 C 3j3 _Ri 

- c 3j4 i?i 4 c 4j3 - 2c 4)4 i?2 4 c 3 , 4 -Ri 4 c 4j3 c 3i3 - c 3>4 i?2 4 c 4 , 3 (25) 
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FIG. 1. a) A map of the Mediterranean Sea. In the small box we show the area under study, 
b) Shows one of the altimetry images obtained from the satellite. 
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FIG. 2. Linear- log plot of the fraction of variance of the eigenvalues in terms of their index. 
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FIG. 3. Temporal modes associated the first two eigenvalues. Circles are the data for the first 

and squares for the second. 
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FIG. 4. Power spectrum of the first four temporal modes. Circles, first mode; squares, second 
mode; triangles, third mode; and diamonds, fourth mode. 
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FIG. 7. Temporal evolution of the third and fourth temporal modes obtained by numerical 
integration Eq. 02 
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a) Dec. 1992 




b)Jan. 1993 




c) Feb. 1993 





e) End of Marc. 1993 
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f) Beginning of Ap. 1993 




g) End of Ap. 1993 



d) Beginning of Marc . 1 993 h) May 1 993 

FIG. 8. Numerical reconstruction of vortex shedding and motion near the Algerian coast. After 
h) the same sequence is repeated. The lighter area (higher altimetric values) next to the coast 
corresponds to an anticyclonic vortex (the moving coherent structure under study). In the figure, 
it always appears with its corresponding cyclonic vortex (darker area). 
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